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Abstract. We investigate the dynamical evolution of a Nas cluster embedded in Ar matrices of various 
sizes from N = 30 to 1048. The system is excited by an intense short laser pulse leading to high ionization 
stages. We analyze the subsequent highly non-linear motion of cluster and Ar environment in terms of 
trajectories, shapes, and energy flow. The most prominent effects are : temporary stabilization of high 
charge states for several ps, sudden stopping of the Coulomb explosion of the embedded Nag clusters 
associated with an extremely fast energy transfer to the Ar matrix, fast distribution of energy throughout 
the Ar layers by a sound wave. Other ionic-atomic transfer and relaxation processes proceed at slower scale 
of few ps. The electron cloud is almost thermally decoupled from ions and thermalizes far beyond the ps 
scale. 

PACS. 36.40.Gk, 36.40.Mr,36.40.Sx, 36.40.Vz, 61.46. Be Atomic and molecular clusters 



1 Introduction 

Metal clusters have many interesting properties which made 
them a much studied species in the past decades. Their 
valence electron cloud provides a finite fermion system 
with remarkable properties as electronic shell effects and 
strong optical absorption in a narrow frequency band, the 
Mie surface-plasmon resonance, for early reviews see [TJ[2J 
[3]. The resonance depends sensitively on the geometry of 
the cluster and thus provides an ideal handle for analyzing 
and for controlling cluster dynamics. This has been much 
exploited for free clusters in all dynamical regimes, for a 
recent overview see [3] . 

The dynamical scenarios become more involved when 
clusters stay in contact with a substrate, either embedded 
inside or deposited on a surface. The variety of combina- 
tions grows huge. This is generally viewed as an advan- 
tage because it allows to tailor desired properties, e.g., 
for designing bio-markers exploiting the prominent opti- 
cal properties of clusters [5l[6] . Corresponding to the huge 
manifold of possibilities, there is a huge body of publica- 
tions. But comparatively little has been done yet in the 
regime of highly non-linear dynamics induced by intense 
laser fields. 

It is the aim of the present contribution to explore, 
from a theoretical perspective, that non-linear regime. As 
test case, we take Na clusters because these are conceptu- 
ally the cleanest metal clusters, least plagued by interfer- 
ence with core electrons and also the simplest for a the- 
oretical description. The clusters are embedded in an Ar 
substrate. This is an insulator and so we simulate a typ- 



ical scenario for a chromophore in an inert environment. 
Moreover, Ar is a very soft material having a weak in- 
terface energy. It exerts only a very small perturbation on 
the cluster which, in turn, maintains basically its structure 
and optical properties. That feature had been exploited in 
the past for studying optical properties of neutral clusters 
which are otherwise hard to handle as free clusters, see 
e.g. [HE]. The weak-perturbative situation will change if 
we consider violent dynamics where cluster electrons and 
ions are driven to heftier encounters with the surrounding 
substrate atoms. The actual test case will be Nag em- 
bedded in a sufficiently large cavity of the Ar^v system. 
We will consider various Ar system sizes N. This serves 
a double purpose : large N serve as an approximation to 
bulk Ar and comparison with smaller N allows to work 
out typical effects of mixed clusters, which are as such an 
interesting species. We excite the embedded cluster by in- 
tense and short laser pulses to high charge states (Na^ + 
up to Q — 4) and follow the subsequent dynamical evolu- 
tion over up to 8 ps. We will analyze the results in terms of 
detailed views, energy transfers, evolution of global shape, 
and relaxation times. 

Clusters in contact with a substrate are much more 
complex systems than free clusters. The expense for a 
fully detailed theoretical description grows huge. Such cal- 
culations are undertaken where details count, mostly for 
studying structure of small compounds (see e.g [9l [T0lfTl~j 
for deposited and embedded metal clusters in various ma- 
terial combinations). But these subtle models are hardly 
applicable to a truly dynamical situation and the enor- 
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mous effort is not really needed for the sort of exploratory 
studies which we have in mind here. The large differ- 
ence between metal clusters (reactive) and rare gas (in- 
ert) suggests a hierarchical model where the environment 
is handled at a lower level of description. Such approaches 
are widely used in theoretical chemistry under the label 
Quantum-Mechanical/Molccular-Mechanical model (QM/MM), 
see e.g. [TJ], and in surface physics, see e.g., [13] ■ A model 
in that spirit for Na clusters in contact with Ar was pro- 
posed in [14] and applied to structure, optical response 
and moderately non-linear dynamics of embedded clus- 
ters in [T5l[T6] . Here we are going to carry forth the stud- 
ies in the regime of strong laser excitations inducing high 
charge states and large rearrangements of cluster and ma- 
trix structure. 



The energy of the Na cluster i/Naciuster consists of TDLDA 
(with a self-interaction correction, see Sec. 12. 2p for the 
electrons, MD for ions, and a coupling of both by soft, 
local pseudo-potentials; that standard treatment is well 
documented at many places, e.g. [4ll8|. The Ar system 
and its coupling to the clusters are described by 



£a, = ]T 



2M a , 



drp A 



(2) 



2 Formal framework 

2.1 Short summary of the model 

In order to allow for a sufficiently large Ar matrix, we use 
a hierarchical approach. We sketch it here briefly and refer 
to [MlTTf] for a detailed layout. 

The Na cluster is treated in full microscopic detail 
using quantum-mechanical single-particle wavefunctions 
{cp n (r, t),n = 1, . . . , N e \} for the valence electrons. These 
are coupled non-adiabatically to classical molecular dy- 
namics (MD) for the Na + ions which are described by 
their positions {R/,7 = 1, ... ,N^ a }. In the present pa- 
per, we have N c \ = N^ a = 8. The electronic wavefunc- 
tions are propagated by time-dependent local-density ap- 
proximation (TDLDA). The electron- ion interaction in the 
cluster is described by soft, local pseudo-potentials. This 
TDLDA-MD has been widely validated for linear and non- 
linear dynamics of free metal clusters [4]fT8] . 

Two classical degrecs-of-frecdom are associated with 
each Ar atom : center-of-mass {R a ,a = 1,...,N}, and 
electrical dipole moment which is parameterized as {R' Q = 
R a + d a , a = 1, . . . , N}. Note that the Ar dipole is prac- 
tically handled by two constituents with opposite charge, 
positive Ar core (at R a ) and negative Ar valence cloud (at 
R a ). With the atomic dipoles, we explicitely treat the dy- 
namical polarizability of the atoms through polarization 
potentials [H5]. Smooth, Gaussian charge distributions are 
used for Ar ionic cores and electron clouds in order to 
regularize the singular dipole interaction. We associate a 
Gaussian charge distribution to both constituents having a 
width of the order of the 3p shell in Ar, similar as was done 
in [20] . The Coulomb field of the (softened) Ar dipoles pro- 
vides the polarization potentials which are the dominant 
agents at long range. The Na + ions of the metal cluster 
have net charge <?Na = +1 and interact with the Ar dipoles 
predominantly by the monopole moment. The small dipole 
polarizability of the Na + core is neglected. The cluster 
electrons do also couple naturally to the Coulomb field 
generated by the atoms. 

The model is fully specified by giving the total energy 
of the system. It is composed as 
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The calibration of the various contributions is taken from 
independent sources, except eventually for a final fine tun- 
ing to the NaAr dimer modifying only the term WeiAr of 
Eq.([5]). The parameters are summarized in table [T] The 
third column of the table indicates the source for the pa- 
rameters. In the following, we report briefly the motiva- 
tions for the choices. 
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fit to NaAr dimer 


T/(core) 
^ArAr 


A Ar =1.367 x 10 -3 Ry , # Ar =6.501 a 


fit to bulk Ar 


"^ArNa 


/3 Na = 1.7624 a c 7 1 , a Na = 1.815 a , ^ N a= 334.85, 
CNa,6= 52.5 ao 6 , CWa,8= 1383 ao 8 


after [21] 



Table 1. Parameters for the various model potentials. See text for detail. 



Most important are the polarization potentials defined 
in Eq.Q. They are described by the model of a valence 
electron cloud oscillating against the rare gas core ion. 
Their parameters are : qAi the effective charge of valence 
cloud, mAr = <7Ar"ici the effective mass of valence cloud, 
&A r the restoring force for dipoles, and OAr the width of 
the core and valence clouds. The QAr and &A r are adjusted 
to reproduce the dynamical polarizability arj(w) of the Ar 
atom at low frequencies, i.e. we choose to reproduce the 
static limit aoif^ — O) and its second derivative <x.'Ij(u) = 0). 
The width OAr is determined consistently such that the 
restoring force from the folded Coulomb force (for small 
displacements) reproduces the spring constant ^Ar- 

The short range repulsion is provided by the various 
core potentials. For the Ar-Ar core interaction, Eq.©, we 
employ a Lennard- Jones type potential with parameters 
such that binding properties of bulk Ar are reproduced. 
The Na-Ar core potential Eq.([7]) is chosen according to 
[21] . Note that the Na-Ar potential from [21] does also 
contain a long range part oc OAr which accounts for the 
dipole polarization-potential. Since we describe that long 
range part explicitely, we have to omit it here. We thus 
choose the form as given in VX rNa - 

The pseudo-potential W e iAr, Eq.©, for the electron- 
Ar core repulsion has been modeled according to the pro- 
posal of [20| . Its parameters determine sensitively the bind- 
ing properties of Na to the Ar atoms. We use them as a 
means for a final fine-tuning of the model. The benchmark 
for adjustment is provided by the Na-Ar dimer. The data 
(dimer binding energy, bond length, excitation energy) are 
taken from [2"2"ll23| . The adjustment shows some freedom 
in the choice of A e \. We exploit that to produce the softest 
reasonable core potential. 

The Van-der-Waals energy, £VdW, Eq.© stems from a 
correlation of the dipole excitation in the Ar atom coupled 
with a dipole excitation in the cluster. We exploit that the 
plasmon frequency WMic is far below the excitations in the 
Ar atom. This simplifies the term to the variance of the 
dipole operator in the cluster, using again the regularized 
dipole operator f Q corresponding to the smoothened Ar 
charge distributions. 



2.2 Validity and limitations 

A few words are in order about the range of validity for 
that hierarchical model. The structure of the coupled sys- 
tems is well described by construction, see the predeces- 
sor in [5D] and the extensive testing in [TJ]. The same 
holds for optical response (see [Foll24j ). However, one has 
to remain aware that the model has limitations with re- 
spect to allowed frequencies and amplitudes. The dynami- 
cal response of the Ar atoms is tuned to frequencies safely 
below the Ar resonance peak, i.e. safely below 15 eV. 
This is well fulfilled in the present calculations where the 
highest relevant frequency is the cluster's Mie plasmon at 
around 3 eV. The amplitudes of the dipole oscillations 
in Ar should also stay below the threshold where free 
electrons are released from the Ar atom into the matrix. 
That could become critical for the violent processes stud- 
ied here. We checked the threshold for electron emission 
by fully quantum mechanical calculations of laser excita- 
tion on Ar atoms and we find a critical field strength of 
order of 1.4 eV/ao. That value holds for constant fields. 
The case is more forgiving if the critical field strength is 
exceeded only for a short time. Anyway, we protocol dur- 
ing our calculation the actual field strength at each Ar 
site. We are on the safe side for the cases with charge 
Q = 3 and lower. The dynamics with charge Q = 4 shows 
occasionally field peaks above the limit, however for very 
short times. 

The quantum mechanical treatment of the Na cluster 
involves two approximations which may limit the quanti- 
tative value of the results. The electron cloud is treated 
with axially averaged potentials and the TDLDA is aug- 
mented by a self-interaction correction (SIC) at the level 
of average-density SIC (ADSIC) [25ll26j . Both approxima- 
tions enhance the barriers for fragmentation of the clus- 
ter, or pieces thereof. The axial approximation is needed 
to enable the large scale calculations performed here. The 
SIC is compulsory for an appropriate description of ioniza- 
tion. We pay the price that our calculations provide rather 
a qualitative picture. The effects shown are certainly rel- 
evant. The thresholds concerning charge state and laser 
field strength are probably not too precise. However we 
also dispose of a full three-dimensional treatment for the 
valence electrons which, moreover, can handle various lev- 
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els of SIC and even exact exchange. We used that for oc- 
casional counterchecks and we never found significant (i.e. 
qualitative) deviations from the approximate treatment. 



2.3 Initialization and propagation 

The test cases in the following studies are Nas as free clus- 
ter and embedded in Ar^r where N = 30, 164, 434, and 
1048. The Ar structures contain a cavity equivalent to 13 
Ar places. This was found to provide sufficient space for 
convenient embedding of Nas- The latter consists of two 
rings carrying 4 Na ions each. The rings have the same 
diameter and are twisted by 45° relative to each other in 
order to minimize Coulomb energy. That structure is the 
same in all cases. The effect of the Ar surroundings on 
the structure is very small. For more details on structure 
and basic optical properties see [15] . All calculations start 
from fully relaxed stationary ground state configurations. 
A laser field is applied as time dependent dipole field with 
frequency wi asor = 2 eV and a temporal sin 2 envelope with 
full width at half maximum (FWHM) of 33 fs. The inten- 
sity (= field strength) is varied to reach the different final 
charge stages Q = 2, 3, and 4. It has typical values in 
the range 1-5 x 10 12 W/cm 2 . Its polarization is along the 
symetry axis (z axis) of the system. 

The numerical treatment uses the standard techniques 
of the time-dependent local-density approximation for the 
cluster electrons and molecular dynamics for Na + ions or 
Ar atoms (TDLDA-MD) as outlined e.g. in [41115] . The 
wavefunctions and fields are represented on a grid in co- 
ordinate space. The quantum mechanical time stepping 
is done by the time-splitting method and the molecular 
dynamics by velocity Verlet. We use absorbing boundary 
conditions to simulate correctly electron escape. 



3 Charge stabilization 

The laser irradiation leads to strong ionization of the clus- 
ter. An interesting issue is then the Coulomb-stability of 
the cluster. This has been much studied in the past under 
the key word "appearance size", the critical size above 
which clusters of a certain charge state become ultima- 
tively stable (i.e. appear in a mass spectrograph), for a 
review see [57] ■ It was shown that the value depends sen- 
sitively on the excitation mechanism, with ns laser pulses 
being less favorable, and fs laser pulses as well as fast ion 
collisions producing significantly lower critical values [28[ 
[29] . For Na clusters, the smallest Na cluster with charge 
state Q = 2 has been observed for Na^ + in experiments 
with laser pulses [3D]. The case of free Na^ + , of inter- 
est here, is at the edge of stability. Higher charge states 
Q > 3, however, lead for sure to immediate Coulomb ex- 
plosion of free clusters. The interesting question is then 
how the surrounding Ar environment modifies the charge 
stability. We will address this question first in terms of 
asymptotic stability by energetic estimates and then an- 
alyze detailed dynamical processes. We will see that the 



Ar matrix manages to prevent immediate Coulomb frag- 
mentation for surprisingly high charge states. The strong 
oscillations of the imprisoned clusters will be analyzed in 
terms of global shape parameters. 

3.1 Stability in principle 

TDLDA-MD simulations for free Na^ + clusters lead to 
immediate Coulomb explosion for charge states Q > 3. 
The case of Nag~ + is at the fringe of stability where small 
changes in the conditions produce large effects in the re- 
sults. Monte-Carlo cooling still produces a stable system. 
The results of dynamical simulations starting from neutral 
Nas and fast ionization depend on the level of approxima- 
tion. The calculations with axial approximation for the 
electrons show huge shape oscillations with size varying 
within a factor of two. Fully three dimensional treatment 
(without SIC and mere TDLDA) starts in the same man- 
ner for the first 5 ps and then "discovers" the Na + emis- 
sion path leading to final Coulomb fragmentation. How- 
ever, the detailed propagation and the breakup channel 
depends critically on the initial configuration. Minimal 
shifts of the ions can lead into totally different directions. 
It is in any case obvious that the cluster is at least not in- 
stantaneously unstable. This observation is corroborated 
by computing (with the 3D code) the Born-Oppenheimer 
surfaces in the space of collective deformation (monopolc, 
quadrupole). We find indeed a local minimum for slightly 
oblate Na^ + surrounded by reaction barriers against fis- 
sion, as can be seen from Fig.Q] The lowest barrier in that 
space is about 1 eV. However, the reaction barrier dis- 
appears in symmetry breaking directions with single-ion 




Fig. 1. Born-Oppenheimer surface for free Nag~ + in the plane 
of collective deformations. The energies are shown versus the 
stretching factor in z and x direction. Three dimensional LDA 
calculations were used. 
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Fig. 2. Binding energies for NaJ + clusters, free (open sym- 
bols) or embedded in Ar434 (filled symbols) . The different sym- 
bols characterize the charge state : Q = <-> circles, Q = 1 <-» 
squares, Q = 2 <-> up-triangles, Q = 3 <-» down-triangles. The 
faint dotted lines connect clusters along same charge states Q. 
The arrows show the energetic path for emission of one Na + 
ion. All clusters have been fully relaxed into their optimal con- 
figuration. 



escape. The Coulomb excitation by immediate ionization 
thus produces a metastable state. It requires some time 
to develop sufficient asymmetry for breakup. The exam- 
ple shows that a case like Na^ + which lies at the limits 
of stability imposes high demands on the theoretical de- 
scription. Nonetheless, one can read the result positively: 
all approaches agree in predicting a long initial period of 
quasi-stability which is well in the range of experimental 
observation, e.g. by pump-and-probe techniques assessing 
the time evolution of cluster shape through the time evo- 
lution of the Mie plasmon resonance [3"Tll3"2"Il3"5] . 

Before proceeding to the dynamical simulations, we 
want to calibrate our expectations by looking at the asymp- 
totic stability. To that end, we compute the ground-state 
binding energies of various Na clusters for several charge 
states. The results for free and embedded clusters are 
shown in Fig. [2j The steady down-slope of the dotted lines 
shows that the binding energy for fixed Q increases in al- 
most constant amounts with the number n of Na + ions. 
The energy difference in vertical direction represent the 
(adiabatic) ionization potentials which naturally increase 
with increasing charge states. The most interesting part 
is to look at the decay paths of the highly charged clus- 
ters, going along emission of a Na + ion (arrows). Down 
slope means energy gain and thus asymptotic instabil- 
ity. A system may be locally stable for quite a while (see 
the above discussion). However at long times, the system 
will find its way to the configuration with lower energy. 
Free Na^ + is indeed found to be asymptotically unsta- 
ble. The energy difference shrinks to a negligible amount 
for embedded Na^ + . This system will then be stable. We 
see generally that embedding has a stabilizing effect for 
charged clusters. For example, the case of Na§ + which 
was clearly explosive for the free cluster is "downgraded" 



to a situation which is comparable to free Nag . Energy 
differences alone however are not fully conclusive for the 
stability times. These depend also on the phase space of 
decay channels and on the reaction barriers. The latter are 
certainly larger for embedded clusters, since a single Na + 
would have to find its way through the closely packed Ar 
medium. 



3.2 Detailed view 

Let us now turn towards the analysis of the dynamical sce- 
nario following the irradiation of an embedded Nas- The 
initial reaction to laser excitation is exactly the same as 
we had already observed previously for the more moder- 
ate cases [16]. There is a fast direct emission of electrons. 
The pattern are the same for free and for embedded clus- 
ters. The Ar environment does not impose any hindrance, 
even for the largest system in our sample carrying 1048 
Ar atoms. The fast ionization is completed at the end of 
the pulse, i.e. at 100 fs, and we are left with a cluster in 
a certain charge state Q depending on the laser intensity. 
The finite net charge of our test system is related to the 
finite size of the setup. In a macroscopically large matrix, 
the electrons would be stuck somewhere within the range 
of their mean free path. They may even drift very slowly 
back towards the now attractive cluster, and eventually 
recombine there. We have simulated that by using reflect- 
ing boundary conditions rather than absorbing ones and 
we find recombination times of the order of several 10 ps. 
Thus the present scenario should provide a pertinent pic- 
ture for several ps, the time window studied here. 

The Coulomb pressure generated by the almost instan- 
taneous ionization drives the cluster ions apart which, in 
turn, stir up the Ar atoms. The dominant effect is expan- 
sion. Thus a visualization in terms of radial coordinates 
is most appropriate. This is done in Fig. [3] showing the 
detailed time evolution for the various charge states and 
system sizes. The cavity is spacious, leaving a large ini- 
tial separation of the Na ions from the first shell of Ar 
atoms. Thus we see for the initial 200 fs a fast expan- 
sion/explosion of the cluster almost as in the free case. 
The expansion is stopped instantaneously if the cluster 
runs against the repulse Ar core. The stopping radius in- 
creases with the charge state Q. After stopping, the cluster 
turns over into oscillations for a while. The Ar atoms take 
up the momentum from the stopped ions. That momen- 
tum propagates like a sound wave through the Ar medium. 
The perturbation is strong enough for Q = 4 to produce 
finally a direct emission of Ar atoms when the wave has 
reached the outmost shell. On the way out, the wave dis- 
tributes the energy very quickly over all shells exciting 
them to more or less hefty fluctuations. The Ar amplitudes 
increase with increasing charge state Q and decreasing sys- 
tem size. The charge dependence is related to the initial 
amount of Coulomb energy oc Q 2 . The size dependence 
comes from the sharing of energy over the Ar degrees of 
freedom. More Ar atoms leave less energy per atom and 
subsequently a smaller amplitude of motion. The Ar ma- 
trix appears surprisingly robust for N = 434 and larger. 
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time evolution of radial coordinates 
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Fig. 3. Time evolution of radial coordinates r n = \J x\ + y% + z% of ions (heavy lines) and atoms (faint lines) for embedded 
Nas after irradiation with a short laser pulse having cji ascl - = 2 eV and sin 2 envelope with FWHM = 33 fs. The intensity is 
varied to reach the different final charge stages Q as indicated. The plotting of the Ar atoms is stopped a bit earlier to display 
more clearly the final status of the Na ions. 



Heating suffices for TV = 30, Q = 3 and TV = 164, Q = 4 to 
induce a steady dissolution of the Ar environment within 
a few ps. The process starts from the outer shells (mind 
that melting starts generally at the surface [31] ) such that 
the Na cluster remains still captured by the Ar shells for 
a long time. The motion of the outmost shell shows also 
interestingly large excursions in other cases (TV = 164 at 
Q = 3 and TV = 434 at Q = 4). The high net charge 
of the total system together with the dipole polarizabil- 
ity of the Ar atoms generates a long-range attractive po- 
tential oc r -4 which allows for these large fluctuations of 
still bound atoms. The stabilization of the cluster appears 
rather limited for the highest charge state Q = 4. At about 
4 ps, the cluster expansion revives. It proceeds, however, 



at a very slow time scale. One could call that a Coulomb 
driven diffusion through the medium. Such a process is not 
visible within the studied time span (up to 10 ps, shown 
are 7 ps for graphical reasons) for charge state Q = 3. It 
cannot be excluded on a longer time span. However, in any 
case up to Q = 4, we see a temporary stabilization inside 
the cavity which lasts sufficiently long for an observation. 

We can also read off from Fig. [3] a few time scales. 
The initial Coulomb explosion lasts only for about 200 
fs. The distribution of energy over the Ar shells is done 
within 0.5-0.7 ps. In fact it propagates with the speed of 
the sound wave which is about 30 ao/ps for Q = 3 and 
even faster for Q = 4. This is close to the velocity of 
sound in pure Ar, as we will see in the next section. The 
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Fig. 4. Trajectories in the a;-z-plane for Q — 4 and N = 434. 



temporary stabilization of the cluster lasts about 4 ps for 
Q = 4, a yet unknown longer time for Q = 3. The thermal 
dissolution of the Ar system takes several ps starting from 
outmost shells. 

The instantaneous stopping, the sound wave and the 
subsequent direct emission of Ar atoms are interesting pro- 
cesses which deserve further inspection. Fig.Q]shows a pro- 
jection of ionic and atomic trajectories onto the x-z plane 
(where z denotes the symmetry axis, and laser polariza- 
tion axis). The ions start with a radial Coulomb explosion, 
are bent over with little energy loss by about 135° and 
hit almost head on the next Ar atom. This then triggers 
the sound wave. The wave is not well visible in coordi- 
nate space. It is basically a momentum wave transmitted 
through small kicks of the atoms. The final kick, however, 
releases the atom in the last shell for a long travel. All 
eight ions perform the initial looping similarly and eight 
Ar atoms are eventually released in that case of Q = 4 
and N = 434. 



3.3 Global shape 

Complementing information can be gained from charac- 
terizing the system through global shape parameters. The 
leading quantity is the overall extension which can be 
quantified by the root-mean-square (r.m.s.) radius r. Next 
important is the quadrupole deformation which we de- 
scribe by the dimcnsionlcss quadrupole moment /?2- Both 
observables are defined in detail as 
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Fig. 5. Time evolution of atomic shape parameters (|12l) . 
r.m.s. radius r and deformation 02, for a pure Ar447 cluster 
and NagAr434 at two charge states. 



These shape parameters are computed for the Na cluster 
as well as for the Ar system. The sum runs up to p = 8 for 
the Na cluster and up to p = N, the number of Ar atoms. 

Fig.[5]shows the time evolution of atomic shape param- 
eters for two charge states in Ar 4 3 4 and compares it with 
the corresponding pure Ar cluster Ar 44 7 which contains 
also the 13 atoms of the cavity. The Ar 44 7 was excited by 
the same laser pulse as for the case Q — 3. The laser has 
a handle on the Ar atoms through their dynamical dipolc 
polarizability. The oscillations have surprisingly much in 
common. Besides the amplitudes, the pattern are similar 
and so are the cycle times, for radial oscillations about 
2-3 ps and for quadrupole oscillations about 6-7 ps. The 
large charge at the center of the system seems to change 
very little on the global properties. The matrix is also in 
that sense inert. The cleanly developed modes allow to 
estimate the sound velocity in these large clusters. The 
radial compression mode corresponds to the longitudinal 
sound mode in bulk material. Its frequency is u> v ih w 1.8 
meV. The momentum of the radial wave is q = tt/R where 
R = 30 ao is the cluster radius. The sound velocity is then 
estimated as w S ound = w V ib/<7 ~ 30 ao/ps which is very 
close to the propagation speed as observed in Fig. [3J The 
quadrupole mode is a surface mode. Its lower frequency of 
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about 0.7 meV yields a velocity of 12 ao/ps corresponding 
to a surface sound wave. 

The amplitudes are, of course, different for the vari- 
ous cases (note the vertical scales in Fig. [SJ. The strong 
laser field has a very small, although visible, effect on the 
pure cluster. The strong reaction for the composite sys- 
tem is due to the chromophore embedded at the center. 
This couples strongly to the laser, becomes ionized, and 
transfers a large portion of its excitation to the Ar sur- 
roundings. The strength of the initial excitation (that is, 
ionization) translates directly into the amplitude of the 
indirectly triggered oscillations. 

Fig. [^summarizes the time evolution of ionic and atomic 
shape parameters for a great variety of systems (left col- 
umn with Nai = 434 and varying charge Q, right column 
with Q = 3 and varying Nai)- The ionic radii show again 
the fast stopping. The stopping radius depends strongly 
on the charge state Q (left lower panel) but seems inde- 
pendent of matrix size (right lower panel). What counts is 
the first shell of Ar atoms whose radius is the same in all 
Ar surroundings [T3] . We also nicely see the steady oscilla- 
tions for Q < 3 and the start of Coulomb diffusion at 4 ps 
for Q = 4. Comparing the ionic r rms for Q = 2 and Q = 3 
(lower left panel), we see that larger Q leads to stronger 
damping of the ionic oscillations because the ions, coming 
closer to the Ar atoms, couple more efficiently to the Ar 
system. The Na + -Ar relaxation time which is very long 
for moderate excitations (see Q = 2 here and [16]) shrinks 
to about 4 ps for Q — 3. The effect of system size on 
the relaxation (right lower panel) is small when compar- 
ing N = 434 with 164. The case of N = 30 is too chaotic 
to be conclusive. 

The effect on the atomic radii (second lower panel in 
Fig. [6]) is related to stopping radius (left) and to system 
size (right). Larger ionization translates to larger ampli- 
tudes up to direct atom emission (Q = 4). Smaller system 
size yields more energy per atom and also increases the 
amplitudes (right panel). Particularly impressive is the 
case of N = 164 with Q = 3. This system shows com- 
paratively huge non-linear fluctuations and is thus at the 
onset of melting. 

The quadrupolc moments (upper two panels) behave 
analogously. The Na ions show a strong initial reaction 
and then develop in the average a large oblate deforma- 
tion (negative The damping again strongly depends 
on charge state Q and less on system size. The final frag- 
mentation for Q = 4 goes clearly to the oblate side. Ac- 
tually the Coulomb force enhances deformation and we 
start from the oblate side. That means geometrically that 
the ions depart preferably orthogonal from the symmetry 
axis (and thus the laser polarization). The trend to the 
oblate side carries over to the Ar atoms in unstable sit- 
uations. Still, the average deformation remains negligible 
in the stable cases, showing again the inertness of the Ar 
matrix. 
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Fig. 7. Time evolution of the bond length fluctuation (|13|) . 
Lower panel : Case of Nas Ar434 at three different charge states. 
Upper panel : Charge state Q = 4 for three different matrix 
sizes. 



3.4 Bond lengths 



The shape parameters (fT2|) provide a global measure. A 
complementing view inside the system is provided by quan- 
tifying the bonding structure. It is related to the distances 
between the Ar atoms. These are fixed in a totally frozen 
configuration. The temporal fluctuations of the distances 
help to conclude on the stability of a structure. A conve- 
nient measure is the average bond length fluctuation 



-'bond 



N(N — 1) 



E 



OX 



.1/2 



(13) 



where (. . .) is a time-average over a large interval. We 
use here 1 ps averaging time. The bond length fluctuation 
is dimcnsionless, rescaled by the average bond distance. 
Thus it measures selectively the state of order. A collective 
deformation (radial expansion, quadrupolc deformation) 
thus does not contribute significantly. 

Fig. [7] shows results for the bond length fluctuations. 
They remain surprisingly small in all cases. At first glance, 
this seems to contradict the huge spatial fluctuations seen 
in Fig. [3] Two points are nevertheless to be noted. First, 
the fluctuations start at the outer shells while the major- 
ity of inner Ar atoms remains still better behaved, and 
second, the bond length fluctuations measure the intrinsic 
order while being insensitive to global motion. The re- 
sult indicates that the Ar system stays surprisingly cold 
in spite of the sometimes hefty reactions. This conclusion 
will be corroborated by an analysis of the internal tem- 
peratures in the following section. 
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Fig. 6. Time evolution of atomic and ionic shape parameters, Eq. (|12[) . r.m.s. radius r and deformation /32, for several test 
cases in comparison. Left panel : Variation of charge state Q for matrix size iVAr = 434. Right panel : Variation of matrix size 
for charge state Q — 3. 



4 Energy transfer 

Less visible but equally important observables are the vari- 
ous energies in the subsystems and the flow between them. 
In fact, energy transport and thermalization are most im- 
portant issues in composite systems, e.g. for analyzing ra- 
diation damage in materials [35,36 . The simplest way to 
characterize the share of energies uses the kinetic energies. 
These are one-body operators and allow an unambiguous 
definition for each subsystem, for the cluster electrons, the 
ions, and the Ar atoms. Fig. [8]shows the ionic and atomic 
kinetic energies for Q = 3 and various matrix sizes. The 
initial phase is detailed in the left panels. One sees the ex- 
plosive growth of ionic energies and the sudden stopping 
with nearly complete energy exchange between ions and 
atoms. About 2/3 of the initial kinetic energy is kicked 



over at once to the first Ar shell. The ions turn to oscilla- 
tions and the damping of the oscillations, seen in Figs. [3] 
and[6l leads to a steady reduction of ionic kinetic energy 
The further evolution of the atomic kinetic energy stays 
below the first peak. Part of the initial energy gain is in- 
vested into potential energy for spatial rearrangement of 
atomic positions and dipoles. Something new happens at 
about 2 ps, particularly for N = 434. The kinetic energy 
increases by about 50%. Part of that energy comes from 
ionic relaxation. Another part has to stem from further 
spatial rearrangement into an energetically more favorable 
configuration. After that, the atomic energy stays almost 
constant. All further energy transfer processes are too slow 
to be resolved within simulation time. 
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Fig. 9. Time evolution of electronic, ionic and atomic kinetic temperatures for NagAr434 at three different charge states. For 
the atoms, we show the temperature from total kinetic energy and from "intrinsic" kinetic energy for comparison (the kinetic 
energy from collective radial motion is subtracted). 
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Fig. 8. Time evolution of ionic and atomic kinetic energies 
for Nag in matrices of different size at charge state Q = 3. The 
left panel zooms on the first few 100 fs. 

However total energies are not the best measure for 
internal excitation. Energy per particle is better suited. 
That leads naturally to the concept of the kinetic tem- 
perature Tkin = 2£kin/(3iV). The definition is a bit more 
involved for the electrons. The kinetic energy as expecta- 
tion value of wavefunctions contains an unavoidable offset 
from quantum uncertainty and from the Pauli principle. 
We subtract that offset, as well as possible collective flow 
contribution, to define an intrinsic electron temperature, 
following the procedure described in [18] . 

The various kinetic temperatures are shown for N = 
434 and all three charge states in Fig. [9] Note that elec- 



trons can use the same scale for the different Q. Their 
intrinsic temperature grows with Q but less than linearly. 
That is due to the strong cooling through ionization, see 
the discussion of Fig. [10] later on. In any case, the elec- 
tronic temperatures arc much higher than any other and 
there is no sign of relaxation towards the other parts, on 
the time span explored here. This shows that from a ther- 
mal point of view electrons and ions/atoms are, to a large 
extent, still decoupled. One can conclude that thermal re- 
laxation times seem to be far beyond our analysis. The 
large electronic temperatures call, in fact, for a description 
beyond pure mean field. Electron-electron collisions play a 
role in that regime. These could be included by switching 
to a semi-classical Vlasov-Uehling-Uhlcnbeck description 
of the electron cloud [171I37] . But long-time semi-classical 
propagation is extremely hard to stabilize and, as we have 
seen above, electrons are almost decoupled from ions, ther- 
mally speaking. Thus we can neglect that effect for the 
present study, at least on the time scales explored here. 

The ionic temperatures show somewhat larger varia- 
tion, still fitting almost on the same scale. The initial peak 
(see lower left panel in Fig. [H]) is not fully shown here. It 
grows with a rate between linear and quadratic in Q. The 
further evolution shows the reverse trend. The kinetic en- 
ergies shrink the faster the higher Q. That is due to the 
increasing coupling to the first Ar shell which accelerates 
the energy exchange to the Ar system. 

Totally different temperature scales are required for 
the Ar matrix. Two effects cooperate here : The ionic ki- 
netic energy grows with Q and the coupling to the Ar 
is strongly enhanced with the increasing ionic stopping 
radius, see lower left panel of Fig. [6] The case of Q — 2 
does not show the abrupt stopping and thus does not have 
the large initial jump in energy. The temperature rather 



F. Fehrer et al: Hindered Coulomb explosion of embedded Na clusters 



11 



> 

</> 
o 

<D 

c 



40 
35 
30 
25 
20 
15 
10 
5 

- 



Na 8 Ar 434 



Ar matrix 
intr. elect. 



ionization 



2 3 4 

charge state Q 
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elect." stands for the intrinsic electronic kinetic energy left in 
the remaining electron cloud, "Ar matrix" stands for the ki- 
netic and potential energy stored in the Ar matrix after relax- 
ation (Na + ions have then only a very small share) . 



estimated that for the present case by integrating the ion- 
ization potentials (IP) along the ionization curve and by 
computing the IP as a function of charge for fixed ionic 
configuration. The results in Fig.QjJlprove that indeed the 
by far dominant portion of energy is carried away with the 
directly emitted electrons. That removes typically 80% of 
the laser energy. The next two large contributors are the 
intrinsic electronic energy and the final energy of the Ar 
atoms (which was, to a large extent, initially ionic kinetic 
energy). The share depends strongly on the charge state. 
Very little is left for the Ar atoms at Q — 2 while they 
gather almost equal share for Q = 4. Here we have to 
keep in mind that the energy balance in Fig. [10] repre- 
sents the status at the end of the simulation, i.e. at 6 
ps. Final full thermalization (if reached before dissolution 
of the system) will give the kinetic energy of electrons, 
ions and atoms, with a thermal share proportional to the 
number of particles in each subsystem. We have order of 8 
electrons and ions but much more atoms. Thus the Ar sys- 
tem will acquire finally almost all residual energy which 
amounts to about 3-6 eV, corresponding to a heating of 
about 50-100 K. 



grows steadily in relation to the equally steady cooling of 
the ionic cloud. The relaxation time for the thermaliza- 
tion is of order of 10 ps. That relaxation time drops to 
2-3 ps for the larger Q with their more intimate Na + -Ar 
coupling. Moreover, most of the energy transfer is concen- 
trated in the stopping and comparatively little is left for 
the subsequent energy exchange. Within the resolution of 
the results, we can state that Na + -Ar thermalization is 
achieved after 4 ps for Q = 3. The energy continues to 
increase for Q = 4. That is due to the Coulomb diffusion 
which converts steadily potential into kinetic energy. 

The results from section [3] show a large amount of col- 
lective motion in radial direction and in quadrupole defor- 
mation. That part of motion is still included in the total 
atomic kinetic energy In order to single out the irreg- 
ular part corresponding to truly intrinsic excitation, we 
have computed the collective kinetic energy contained in 
the radial oscillations (the slower quadrupole part con- 
tributes very little) and subtracted it to obtain an intrin- 
sic kinetic energy. The corresponding kinetic temperature 
is also shown as the curves labeled "atoms intrinsic" in 
the lower panels of Fig. [9] The collective part is always 
the largest portion and even dominates by far for Q > 3. 
The minor part is left for intrinsic motion. That explains 
the findings from the bond lengths in Fig. [7] which show 
very little perturbation of the internal structure of the Ar 
system. In spite of the huge amounts of energy flowing 
around, that structure persists astonishingly robust. 

Finally Fig. [10] shows the distribution of the excita- 
tion energy initially induced by the laser. The total height 
stands for the excitation energy. Of course, it increases 
with increasing intensity and subsequently increasing charge 
state. This net investment grows a bit faster than linear. 
Earlier experience with free clusters shows that a major 
cooling mechanism is electron emission [38 ] l39 j . We have 



5 Conclusions 

We have analyzed the laser-induced dynamics of a Nas 
cluster embedded in an Ar matrix. We have used a hier- 
archical model with a detailed quantum-mechanical treat- 
ment of the clusters electron cloud and a classical descrip- 
tion of the Na + ions as well as of the surrounding Ar atoms 
whereby we include the Ar dipole moments to account 
for the dynamical polarizabilitics. The Ar matrix is ap- 
proximated by a large Ar cluster covering several shells of 
atoms. To distinguish finite size effects, wc have used dif- 
ferent sizes from N = 30 to 1048 Ar atoms. The test cases 
provide several interesting aspects : They are by construc- 
tion a model for a cluster embedded in a rare gas matrix, 
they display also interesting features specific for a finite 
composite, and they are prototypes of a chromophore in 
an else- wise inert material. For the excitation mechanism, 
we consider irradiation by a short 33 fs laser pulse with 
intensity of about 10 12 W/cm 2 . 

It was found in earlier studies of structure and mild 
excitations that the Ar matrix induces only very small 
perturbations of the cluster properties. That remains true 
for the electronic dynamics in the more violent regime. 
However, the now released large amplitude motion of the 
cluster ions interferes with the atomic cage. This produces 
dramatic differences of the cluster dynamics as compared 
with the free case. The (meta-)stability of charged clus- 
ters is much enhanced; even a Na g + cluster remains bound 
for about 4 ps, and lower charge states much longer. The 
Coulomb explosion of the highly ionized cluster is abruptly 
stopped by the first shell of Ar atoms and most of ionic 
kinetic energy is transferred to the Ar shell within a few 
10 fs. The transferred energy is then spread very quickly 
by a sound wave all over the Ar matrix. Subsequent re- 
laxation and rearrangement processes between ions and 
atoms proceed at a time scale of several ps. The Coulomb 
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instability takes over within that time span for Q = 4 
and drives a slow Coulomb-driven diffusion of Na + ions 
through the Ar medium. All processes related to shape 
evolution within that time scale are in reach of experi- 
mental observation by pump-and-probc analysis using the 
cluster surface plasmon as handle. 

We have also analyzed the energy balance. The time 
scales for energy transfer confirm the findings for the ionic- 
atomic relaxation as deduced from shape analysis. The 
electron cloud adds by far the most efficient energy loss 
through direct ionization during laser impact. And even 
after large energy outflow with the emitted electrons, a 
large part of excitation energy remains kept in internal en- 
ergy of the electrons because the thermal (i.e. non-collective) 
coupling of electrons to ions and atoms is extremely weak. 
The small remaining fraction of excitation energy goes 
through Coulomb energy (charging) , kinetic energy of Na + 
ions (Coulomb explosion), quickly to the Ar system. In 
spite of the small amount of initial laser energy finally 
transmitted to the Ar matrix, the Na cluster as chro- 
mophore pumps orders of magnitude more energy into the 
system as a pure Ar cluster would be able to extract from 
the same laser pulse. 
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